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Hard spheres interacting through a square- well potential were simulated using two different meth- 
ods: Brownian Cluster Dynamics (BCD) and Event Driven Brownian Dynamics (EDBD). The struc- 
ture of the equilibrium states obtained by both methods were compared and found to be almost the 
identical. Self diffusion coefficients (D) were determined as a function of the interaction strength. 
The same values were found using BCD or EDBD. Contrary the EDBD, BCD allows one to study 
the effect of bond rigidity and hydrodynamic interaction within the clusters. When the bonds are 
flexible the effect of attraction on D is relatively weak compared to systems with rigid bonds. D 
increases first with increasing attraction strength, and then decreases for stronger interaction. Intro- 
ducing intra-cluster hydrodynamic interaction weakly increases D for a given interaction strength. 
Introducing bond rigidity causes a strong decrease of D which no longer shows a maximum as 
function of the attraction strength. 

PACS numbers: 05.10.Ln, 82.70.Dd, 82.70.Gg 



I. INTRODUCTION 

Suspensions of particles with attraction exhibit equi- 
librium states as well as non-equilibrium states like gels 
or glasses that evolve very slowlji^i^ii. Numerical sim- 
ulations were found to be useful for the understanding 
of the structure and the dynamics of such systems^i^. 
The advantage of simulations is that large scale phe- 
nomena may be related to the microscopic trajectories 
of the particles. In recent years computer simulations 
have yielded valuable insight not only into equilibrium 
properties such as cluster size distributions and struc- 
ture factors, but also into the evolution of the system 
during phase separation's^ and gelation^iiSiiiiii^iH. 

Often Monte Carlo methods are used to study struc- 
tural properties at equilibrium and molecular dynamics 
to study dynamics and kinctics^"'-^^. Montc-Carlo meth- 
ods allow one to study relatively large systems and gen- 
erally require less computer time to obtain equilibrium 
states. A main draw back of classical Monte Carlo meth- 
ods is that the definition of time is usually unphysical 
so that the evolution of the system towards equilibrium 
cannot be compared to that of real systems. Molecular 
dynamics simulate the particle displacement more realis- 
tically, but the system size and time scales that can be 
simulated with the current generation of computers are 
still relatively small. 

The simplest model of interacting fluids is an ensemble 
of hard spheres that interact through a square well poten- 
tial. Here we compare two different methods to simulate 
this model system: Brownian Cluster Dynamics (BCD) 
and Event Driven Brownian Dynamics (EDBD)^^. With 
BCD clusters are constructed by forming bonds between 



spheres within each others interaction range with a given 
probability. With this method it is possible to account 
for hydrodynamic interaction within the clusters, though 
not between the clusters. It is also possible to study the 
influence of bond rigidity on the dynamics. This is im- 
portant because in real systems bonds may be more or 
less rigid. With EDBD hydrodynamic interaction is ig- 
nored and the bonds are inherently completely flexible. 
We will show here that for reversibly aggregating sys- 
tems bond rigidity has no influence on the structure of 
the steady state, but has a huge effect on the dynamics. 
In the following we will first describe the two simula- 
tion methods. Then we compare the structure factors 
and the cluster size distributions of homogenous equilib- 
rium states. We will show that almost the same struc- 
tures are obtained at steady state with both methods. 
The main part of paper deals with the self diffusion coef- 
ficient as a function of the interaction strength. We com- 
pare the results obtained by the two simulation methods 
and discuss the influence of bond flexibility and intra- 
cluster hydrodynamic interaction. 



II. SIMULATION METHOD 

We simulate hard spheres interacting through a square 
well potential characterized by a well depth u and a well 
width e using BCD and EDBD. Both simulation meth- 
ods start with an ensemble of Ntot randomly distributed 
spheres with diameter equal to unity in a box of size L so 
that the volume fraction is defined a,s (f) — {it / 6) Ntot/ L^- 

Event Driven Brownian Dynamics. This method 
was described in detail in the literatures^ and we only 



resume here the principal features. Initially a random 
velocity is assigned to each sphere from a Gaussian dis- 
tribution with average squared velocity < v"^ >= 3kT/M 
and variance {kT/MY^^, where k is Boltzmann's con- 
stant, T is the absolute temperature and M is the mass 
of the particle. Events are defined as occurrences when 
the sphere is at a distance unity or 1 -I- e from another 
sphere, i.e. when spheres touch, or cross the interaction 
range. First the shortest time. At, before an event occurs 
is calculated. Then all spheres are moved over a distance 
r — wAt, where At cannot exceed a certain maximum 
value. It is important that the maximum value of At is 
sufhciently small so that the motion is Brownian over the 
relevant length scales, i.e. the interaction range and the 
average distance between the nearest neighbors. 

After each sphere has been displaced the velocity of the 
spheres involved in the event are changed while conserv- 
ing the energy and the momentum. When the event is 
a collision the spheres will bounce elastically in opposite 
directions. When the sphere enters a well its velocity 
is decreased because the potential energy is increased. 
When the sphere tries to leave a well it either bounces 
back elastically or it exits with a higher velocity. The 
change in the velocity and the probability to exit the 
well depend on u. 

The mean squared displacement of a single sphere is 
given by: < r^ >= n(At)^ < v^ >, where n is the num- 
ber of simulation steps and At is equal to the maximum 
time step. If time is defined as i = n/(Ai)^ then the 
free diffusion coefficient of a single sphere is equal to: 
Do ~ kT/{2M). In this paper the unit of energy is the 
thermal energy. We note that in the literature often u is 
fixed at unity and T is varied. 

Brownian Cluster Dynamics. Spheres are consid- 
ered to be in contact when they are within each others 
interaction range, i.e. when the center to center distance 
is smaller than 1 -f e. In the so-called cluster formation 
step, spheres in contact are bound with probability P. 
Alternatively, bonds are formed with probability a and 
broken with probability /3, so that the P = a/ {a + fi). 
In the latter case one can vary the kinetics of the aggre- 
gation from diffusion limited (a = 1) to reaction limited 
(a -^ 0) with the same P and thus the same degree of 
reversibility. Clusters are defined as collections of bound 
spheres and monomers are clusters of size one. After this 
procedure Nc clusters are formed. We mention that more 
complex interaction potentials can be simulated by mak- 
ing P a function of the distance between two spheres. 

The ratio of the number of bound [vf,) to free con- 
tacts {v-viy) is given by the Boltzmann distribution: 
Vb/{v — Vb) = exp(AiJ), where AH is the enthalpy dif- 
ference between bound and free contacts. The forma- 
tion of Vb randomly distributed bonds over v contacts 
leads to a decrease of the free energy equal to u per con- 
tact. This decrease may be written as the sum of the 
decrease of the enthalpy and the gain of the entropy: 
V ■ u — VbAH — TAS. The latter is determined by the 
number of ways v^ bonds can distributed over v contacts: 



TAS = \a{v\l[vi,\{v - i/fe)!]). Noticing that V^Vb/v, we 
can express P in terms of u: 



P = I — exp(u) 



(1) 



The cluster construction step is followed by one of three 
different movement steps that each simulates a different 
type of cluster dynamics. 

BCDl. Ntot times a sphere is randomly selected and 
an attempt is made to move it a distance s in a random 
direction. The movement is accepted if it does not lead 
to overlap with any other sphere in the system and if it 
does not lead to the separation of bound spheres beyond 
the interaction range. Again it is important to choose the 
step size s sufficiently small so that the motion is Brow- 
nian over the relevant length scales. We have found that 
the results on the equilibrium structure were independent 
of the step size if s was at least five times smaller than the 
interaction range and at least three times smaller than 
the average distance between nearest neighbors^iii. 

The mean squared displacement of a single sphere is 
given by: < r^ >= ns^ where n is again the number of 
simulation steps. Time was defined as t = n/s^, so that 
the free diffusion coefficient of single spheres is equal to: 
Do = 1/6. 

BCD2 Nc times a cluster is randomly selected. An 
attempt is made to move the whole cluster over a distance 
s/d in a random direction with d the diameter of the 
cluster. By definition this cooperative movement never 
leads to bond breaking. The movement is refused if it 
leads to overlap of any of the spheres in the clusters with 
other spheres in the system. The free diffusion coefficient 
of single spheres is thus still 1/6, but the free diffusion 
coefficient of clusters is l/(6d). 

BCD3 This movement step is a combination of the 
previous two. First the movement step BCDI is executed 
and the displacement of the centers of mass of each clus- 
ter is calculated. Then each cluster is given an additional 
displacement in the same direction so that the total dis- 
placement of the center of mass is the same as would 
be obtained by the movement step of BCD2. Again dis- 
placements are refused if they lead to overlap. As for 
movement step BCD2, the free diffusion coefficient of sin- 
gle spheres is 1/6, and that of larger clusters is \/{&d). A 
lower degree of flexibility can be simulated by perform- 
ing movement step BCDl with a lower frequency than 
movement step BCD2. 

The methods EDBD, BCDl and BCD3 simulate sys- 
tems with flexible bonds, while BCD 2 simulates systems 
with rigid bonds. Using EDBD and BCDl the effective 
friction coefficient of clusters is proportional to their ag- 
gregation number (so-called Rouse dynamics) ^^, while 
for BCD2 and BCD3 it is proportional to their radius 
(so-called Zimm dynamics)^^. It is, of course, straight- 
forward to modify BCD2 to simulate systems with rigid 
bonds in which the friction coefficient is proportional to 
their aggregation number. This has not been done here, 
because in reality hydrodynamic interaction causes the 



free diffusion coefficient of clusters in solution to be in- 
versely proportional to their radius??. The movement 
steps of EDBD and BCDl are similar and one expects 
that diffusion coefficients obtained by these methods are 
close. 

BCD may be considered a Monte Carlo type simula- 
tion, but if one is interested only in the structural proper- 
ties at equilibrium it is more efficient to use other Monte- 
Carlo techniques that do not yield realistic kinetics or 
dynamics. BCD does not fulfill the condition of detailed 
balance, but does lead to a steady state independent of 
the starting configuration, which shows that it fulfills the 
condition of balance^i. The same steady state is reached 
for each of the three movement steps with one exception: 
BCD2 does not lead to crystallization. The reason is 
that the pathway to form crystals is extremely unlikely 
when the bonds are rigid. BCD2 is therefore an excel- 
lent method to explore to properties of attractive spheres 
while avoiding crystallization^. 



1.2 



1.0 



0.6 



0.4 



0.2 



0.0 




0.01 



0.1 



III. RESULTS AND DISCUSSION 
A. Equilibrium structure 

The strength of the attraction is determined by the 
well width and the well depth. However, the equilibrium 
structure obtained at different e is close if u is chosen such 
that the second virial coefficient is the same, especially 
if e is small^i^^. B2 may be written as the sum of a 
repulsive (excluded volume) part and an attractive part: 

B2 = Bre 



Batt, with 

B„tt 



P 



[l-P] 



((1+ef-l) 



where the unit of B2 is the volume of a sphere. 

As mentioned in the introduction strong attraction be- 
tween the spheres leads to phase separation, while weak 
attraction leads to a homogeneous equilibrium state con- 
taining a distribution of transient clusters. We have char- 
acterized these states in terms of the static structure fac- 
tor {S{q)) and the cluster size distribution. 

The structure factor at 5 —> is inversely propor- 
tional to the compressibility and can be expressed in 
terms of a virial expansion at small volume fractions: 
1/S'(0) = 1-^2^2^ + 35302. For hard spheres interacting 
with a square well potential the second and third virial 
coefficients have been calculated analytically^^,. Fig- 
ure([T]) shows a comparison of 1/5(0) as a function of 
(j) obtained from BCD simulations with the values calcu- 
lated using the virial expansion. There is good agreement 
up to (/) = 0.1 beyond which higher order virial terms be- 
come important. 

Figure(l2]) compares S{q) obtained with BCD and 
EDBD at several conditions. Within the uncertainty 
range the same structures were observed with the two 
methods. 



FIG. 1: Concentration dependence of the compressibility for 
equilibrium systems at e = 0.5 for B2 = —6 obtained by BCD. 
The solid line represents a calculation using the second and 
third virial coefficients, see text. The error bars represent the 
95% confidence based on the results of 8 simulations. 
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FIG. 2: Comparison of the structure factor of equilibrium 
systems at = 0.15 and e — 0.1 obtained by BCD (filled 
symbols) and EDBD (open symbols) with B2 = 2 (squares) 
or B2 = —2 (circles). 



The cluster size distribution represents a more pre- 
cise characterization. Clusters of bound spheres can be 
formed by connecting spheres in contact with probabil- 
ity P = 1 — exp{u) as defined above. A detailed analysis 
of these cluster size distributions has been reported else- 
where. Here, we have analysed the size distribution of 
clusters formed by connecting all contacts in order to fa- 
cilitate comparison between BCD and EDBD. Figure® 



shows a comparison of N{m) i.e. the average density of 
clusters consisting of m spheres at </> = 0.15 and e = 0.1 
for two values of i?2: 2 and 0. The width of the distribu- 
tion increases with increasing u until at a critical value a 
system spanning transient network of spheres in contact 
is formed. There is a very small, but systematic differ- 
ence between the cluster size distributions obtained with 
BCD and EDBD. Slightly larger clusters are formed with 
BCD. Nevertheless, we may conclude from these exam- 
ples and similar comparisons at other conditions that the 
equilibrium structures obtained by BCD and EDBD are 
almost the same. 



glass transition, which has been the subject of intensive 
investigatioii22i2I. 




FIG. 3: Comparison of the cluster size distribution of equilib- 
rium systems at (p = 0.15 and e = 0.1 obtained by BCD (filled 
symbols) and EDBD (open symbols) with B2 = 2 (squares) 
or 52 = (circles). 



IV. DIFFUSION 

It can be shown that random displacement with con- 
stant step size as done with BCD leads to Brownian dif- 
fusion at distances much larger than the step siz o^^'^^ . 
For a proper comparison of the dynamic properties of 
BCD and EDBD simulations one has to ensure that the 
time scales and the free diffusion coefficients are the same, 
which can be done by choosing At = s and kT/M =1/3. 
Figure ([3]) shows a comparison of the average mean square 
displacement of hard spheres obtained from BCD and 
EDBD at 0=0.3. In both simulations <r2> becomes 
proportional to t and the diffusion coefficient can be cal- 
culated as D =< r^ > /6t. Figure([5]) shows that the 
same (f> dependence of D is obtained by the two methods 
within the uncertainty of the simulations. D decreases 
with increasing volume fraction due to steric hindrance 
and the diffusion slows down critically at the so-called 




FIG. 4: Comparison of the mean square displacement of non- 
interacting hard spheres obtained by BCD (filled symbols) 
and EDBD (open symbols). 

When we introduce attraction between the spheres we 
need to consider cooperative cluster motion and bond 
flexibility. EDBD allows only one type of motion, but 
using BCD one can choose between different movement 
steps. 

We have simulated the mean square displacement of 
spheres using EDBD and BCD with the three different 
movement steps described above. In each case diffusive 
motion was observed for large t and the diffusion coeffi- 
cient could be determined. BCDl and EDBD simulate 
the same situation and therefore the results should be the 
same. Figure® shows an example of the dependence of 
D on the step size obtained using BCDl. It appears that 
the exact value of D is more sensitive to the step size 
than the static structure factor or the cluster size distri- 
bution since the latter did not depend significantly on the 
step size for s < 0.5. The value extrapolated to s = 
is the same as the value found with EDBD within the 
simulation error. Similar results were obtained at dif- 
ferent volume fractions and interaction strengths. The 
fact that these very different simulation methods lead to 
the same results, strengthens confidence in both meth- 
ods. In terms of computational efficiency both methods 
are equivalent. 

A comparison of the dependence of D on Batt obtained 
with BCD using the 3 different movement steps is shown 
in fig.© for two different volume fractions (0.15 and 
0.49) and two different well widths (0.1 and 0.5). The 
range of Batt that can be explored is limited by the liquid- 
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FIG. 5: Comparison of the decrease of the diffusion coeffi- 
cient with increasing volume fraction for non-interacting hard 
spheres obtained with BCD (filled symbols) and EDBD (open 
symbols}^. 



liquid or liquid-crystal phase separation that occurs at 
strong attraction. The values of D shown in fig.® were 
obtained at a single value of s, but for a few examples 
the effect of s was determined, which showed that they 
were about 10% smaller than the values extrapolated to 
s = 0. 

D decreases strongly with increasing attraction when 
the bonds are rigid (BCD2). In this case the displace- 
ment of bound spheres is equal to that of the center of 
mass of the clusters to which they belong. The size of 
the clusters increases rapidly with increasing attraction 
and beyond a critical value a transient (bond) percolating 
network is formed. Spheres that are part of the network 
are immobile until the bonds that tie them to the network 
are broken. A detailed study of the diffusion coefficient 
of square well fluids forming rigid bonds using BCD2 has 
been reported recently'^". It was shown that D decreases 
with increasing Batt following a power law for large Batt 
and only becomes zero when the bonds are irreversible, 
i.e. Batt -^ oo. 

In comparison, the influence of attraction on D is weak 
when the bonds are flexible, i.e. using EDBD, BCDl or 
BCD3. The difference between methods BCDl (EDBD) 
and BCD3 is that for the latter clusters move faster 
(Zimm dynamics) so that D is slightly larger. The effect 
increases with increasing cluster size and is expected to 
be maximal close to the percolation threshold. The val- 
ues of Batt at the bond percolation thresholds are about 
6 at = 0.15 for both well widths, while at (f) — 0.49 
they are 0.5 and 1.2 for e = 0.1 and 0.5, respectively. 
The difference between the two methods decreases for 
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FIG. 6: Dependence of the diffusion coefficient on the step 
size simulated using BCDl (</) = 0.49, Batt = 6, e = 0.1). The 
dashed line represents the result from EDBD. We note that at 
(/!> = 0.49 the average distance between randomly distributed 
spheres is 0.0142^12^. 



larger Batt when most spheres belong to the percolating 
network that has no center of mass movement. The few 
remaining free spheres are mostly monomers so that the 
movement steps BCDl and BCD3 become equivalent. 

Regardless of the method, the effect of attraction on D 
is qualitatively different if the bonds are flexible, because 
in that case bound spheres can move freely within the 
interaction range. D increases weakly with increasing 
attraction until it reaches a maximum beyond which it 
decreases. The relative amplitude of the increase is very 
small for the volume fractions studied here, but becomes 
important at higher volume fractions^^>. It is at the origin 
of the so-called reentrant glass transition of interacting 
spheres as a function of the attraction strengtb^i^. The 
influence of attraction on the critical slowing down of 
hard spheres has attracted a lot of attention in the recent 
past and has been investigated for square well fluids using 
EDBD simulations^. 

The appearance of a maximum diffusion coefficient can 
be qualitatively understood by considering two opposing 
effects. On one hand, attraction causes clustering of par- 
ticles so that more space is created in which monomers 
and clusters can diffuse freely, leading to faster diffusion 
of the spheres. On the other hand, bonds restrict the 
motion of spheres and the long time diffusion of bound 
spheres is equal to the center of mass diffusion of the 
clusters to which they belong. The restriction becomes 
more important as the average bond life time increases. 




FIG. 7: Variation of the diffusion coefficient with increasing 
attraction obtained from BCDf (circles), BCD2 (triangles), 
and BCD3 (squares) at two different volume fractions and 
two different interaction widths. The interaction strength is 
expressed in terms of the attractive part of the second virial 
coefficient. 



When the attraction is weak the average bond life time 
is still small so that the effect of restriction is weak and 



the effect of creating more free space dominates leading 
to an increase of D. With increasing Batt the clusters 
become larger and the average bond life time increases 
until the effect of increasing restriction of the movement 
becomes more important than the effect of increasing free 
volume so that D decreases. These features are indepen- 
dent of the volume fraction and the well width. The effect 
of attraction on the diffusion coefhcient remains small in 
the single phase regime if the bonds are flexible atleast 
for 6 < 0.5. 



V. SUMMARY 

BCD and EDBD simulations of hard spheres interact- 
ing with a square well potential lead to steady states that 
have almost the same structure factor and the cluster size 
distribution. 

EDBD assumes flexible bonds and ignores hydrody- 
namic interaction. The values of the self diffusion coeffi- 
cient obtained by EDBD are very close to those obtained 
with BCD if the same assumptions are used. A weak 
maximum of D is found as a function of the interaction 
strength caused by the opposing effects of increasing free 
volume and increasing bond life time. 

The effect of intra-cluster hydrodynamics (Zimm dy- 
namics) and bond rigidity can be explored with BCD. 
Introducing rigid bonds leads to a strong decrease of D 
with increasing attraction and suppresses the maximum. 
Introducing intra-cluster hydrodynamics to the system 
with flexible bonds weakly increases I? at a given inter- 
action strength. 
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